! 2010 document
! This routine controls the calculation of all canopy gas-exchange.
! It calls the routines GEOM, BALENG and REPORT.
! LAI is here calculated.


!cgfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!c     gfx.for                      V4.0.4                      02.12.97  gy
!c==========================================================================
!c     original Version: Fortran-Version 3.0.7 = C-Version 3.0.7-0.1
!c                       with Ron Ryel's work who included penumbral
!c                       effects at BITOEK in Nov.1995.
!c==========================================================================
!c     V4.0.0: Penumbral effects can be choosen by setting the common-
!c               -variable penu to 0 or 1.
!c               penu=0: no penu. effects       penu=1: with penu. effects
!c             Analytical version can be choosen by setting the common-
!c               -variable anal to 0 or 1.
!c               anal=0: iterative solution     anal=1: analytical solution
!c             Included Eva's p_flag (area basis of rates) and a_flag
!c               (leafarea basis). See comments in subroutine psn6.
!c             avggl(nlay) in subroutine report calculated from gl(nlay,sun)
!c               and gl(nlay,shd), not from transpiration.
!c     V4.0.1: Eva: (cint1-gamma) in psn6 changed to cint1
!c     V4.0.2: gy: data for penumbral is read in from parameterfile
!c               penu=0: no penu. effects       penu>0: with penu. effects
!c               penu is number of elements of the 2nd dimension of yint
!c               and slope in commonblock penudata (z.B. slope(6,penu))
!c               tested with beech and grasland
!c               Eva: 1.56 in psnx changed to 1.6
!c               in guessgpress: itlimit=lower limit of iteration
!c                               (gmin earlier)
!c               fac in psn5 and psn6 is read in from parameterfile as rdfac
!c                   (if rd-parameters are from meassured data during night
!c                    rdfac=0.5 else rdfac=1.0)
!c     V4.0.3: subroutine sunsky: W.Eugster's strst calculation and Ron's
!c             changes.
!c             Bug fixed in rdSeason if ndiam>0 (12,13)
!c     V4.0.4: a_flag, p_flag stuff corrected, rdfac in inputfile moved to
!c             physiology parameters
!c     F90 V0.9 By Bumsuk

!c==========================================================================
      subroutine initSeason
!c
!c     read site specific seasonal constants
!c--------------------------------------------------------------------------
      IMPLICIT NONE

      character*15 filename

      filename='RICE_odae.gfx'

      open(unit=95,file=filename,status='old')
      call rdSeason(95,1)
      close (unit=95)

      return
      end




!cgfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine rdSeason(fp,veg)

!c     read site specific seasonal constants
!c--------------------------------------------------------------------------
      implicit none

      integer drow,dcol,dveg

      parameter (drow=6)
      parameter (dcol=18)
      parameter (dveg=1)

      common/rdsea/matrix(dcol,drow,dveg)

      integer fp,ivar(1),ierr,xneedle,veg
      double precision matrix

!c     Veg.typ (NEEDLES)  A_FLAG  P_FLAG  maxLAI  maxSAI OMEGA
      call read0(fp,matrix(1,1,veg),ivar,6,2,ierr,*99)

!c     NEEDLES
      xneedle=matrix(1,1,veg)

!c     FVC       C      D     O2    ALPHA   RDFAC
      call read0(fp,matrix(1,2,veg),ivar,6,2,ierr,*99)

!c     HA HD DS E EATAU FTAU EAKO FKO EAKC FKC EAVC HDVC DSVC GMIN GMAX GFAC
      call read0(fp,matrix(1,3,veg),ivar,16,2,ierr,*99)

!c     Leaf Angle (degrees), Portion of leaves that is alive, Avg. width of leaves (cm), Stem Inclination Angle
      call read0(fp,matrix(1,4,veg),ivar,4,2,ierr,*99)

!c     bunching factor only for needles
      if (xneedle.eq.1) call read0(fp,matrix(1,5,veg),ivar,2,2,ierr,*99)

      return

 99   if (ierr .lt. 0) then
         print*,' PARAMETER-FILE INCOMPLEATE OR MISSING '
         print*,' UNIT : ',fp
      elseif ((ierr.gt.0).and.(ierr.lt.200)) then
         print*,' ERROR WHILE READING FROM PARAMETER-FILE IOS=',ierr
         print*,' UNIT : ',fp
      else
         print*,' CORRUPT OR INCOMPLEATE LINE IN THE PARAMETER-FILE '
         print*,' UNIT : ',fp
      endif
      close (unit=fp)
      stop

      end





!cgfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       subroutine read0(lunr,fvar,ivar,m,ir,ierr,*)
!c
!c      c:  R. Siegwolf
!c--------------------------------------------------------------------------
!c      THIS SUBROUTINE PROVIDES A PROPER READING OF INDEXED VARIABLES.
!c      IF AN ERROR OCCURES THE READ STATEMENT DOESNT JUMP TO THE NEXT
!c      RECORD AND THUS AVOIDES AN UNWANTED FORWARD POSITIONING IN THE
!c      RECORD OF THE READ FILE.
!c      THIS WAS NECESSARY, BECAUSE WITH APOLLO - AND MICROSOFT FORTRAN
!c      THE CONUNTER ON AN INDEXED READ STATEMENT  IS REINITIALIZED TO
!c      THE INITIAL VALUE AFTER ENCOUNTERING AN ERROR, WHICH IS NOT
!c      THA CASE WITH VAX-FORTRAN AND SOME IBM COMPUTERS
!c      LUNR .... UNIT NUMBER IN THE READ STATEMENT
!c      FVAR .... ARRAY FOR REAL VARIABLES
!c      IVAR .... ARRAY FOR INTEGER VARIABLES
!c      M .... NUMBER OF VARIABLES TO BE READ
!c      IR .... INDICATING WEATHER AN INTEGER (IR=1) OR A REAL (IR=2)
!c              SHOULD BE READ
!c      IERR .... IOSTAT VALUE IF A READ-ERROR OCCURES. IF THE ERROR OCCURES
!c          BY READING FROM THE FILE THE IOS VALUES ARE UNCHANGED IF THE ERROR
!c          OCCURES BY READING THE STRING A VALUE OF 1000 IS ADDED TO IOS
!c      * .... ALTERNATIVE RETURN. A LABEL MUST BE SPECIFIED IN THE CALLING
!c             ROUTINE WHERE THE PROGRAMM MUST RESUME EXECUTION AFTER
!c             ENCOUNTERING A READ-ERROR
!c--------------------------------------------------------------------------
      implicit none

      dimension fvar(*),ivar(*)
      character *512 istr,cn*20,ct*1,fmtf*8,fmti*5
      double precision fvar
      integer ivar,m,lunr,ir,ierr,k,i1,iv,kip,i0,ip,ix,ios,ivnr

      data fmtf /'(f14.14)'/
      data fmti /'(i17)'/
      ct=' '
      cn='go                  '
      k=0
      i1=1
      iv=0
      kip=0

!c     Reading a line form any file as a string ISTR
 1    istr(1:)=' '
      read(lunr,'(a)',iostat=ios,err=99,end=99)istr
      ix=index(istr(1:5),'#')
      if (ix.ne.0) then
         goto 1
      endif

!c     The string ISTR is parsed for integer and real variables and
!c     the character sequence is stored in the substring CN. After
!c     detecting a blank or comma as separator the format variables
!c     fmtf = for real and
!c     fmti = for integer variables is adapted to
!c     the lenght of the variable. Then the substring is converted to
!c     a real or integer variable and stored into the corresponding
!c     array, which can be returned to the calling procedure.
!c     ip = indicates the position of the decimal point in the
!c     substring cn.
!c     kip = is the format descriptor after the decimal point
!c     k = number of characters for the variable
!c     ct = character variable for parsing the string ISTR
!c     i1 = number of variables decoded
!c     iv = is set to 1 if a character other then blank or comma
!c     is encountered and is reset to 0 if comma or blank is
!c     found

      do 2 i0=1,512
!c        check for empty lines in data file
         if ((i0.gt.500).and.(i1.eq.1)) goto1
         if (i1.le.m) then
            read(istr(i0:i0),'(a1)') ct
            if ((iv.eq.0).and.(ct.eq.' '.or.ct.eq.',')) goto 2
            if (ct.ne.' '.and.ct.ne.',') then
               k=k+1
               cn(k:k)=ct
               iv=1
            else
               if (ir.eq.1) then
                  write(fmti(3:4),'(i2)') k
                  read(cn(1:k),fmti,iostat=ios,err=90) ivar(i1)
               endif
               if (ir.eq.2) then
                  ip=index(cn,'.')
                  if (ip.gt.0) kip=k-ip
                  write(fmtf(3:4),'(i2)') k
                  write(fmtf(6:7),'(i2)') kip
                  read(cn(1:k),fmtf,iostat=ios,err=90) fvar(i1)
               endif
               i1=i1+1
               cn='go                  '
               k=0
               iv=0
               kip=0
               goto 2
            endif
         endif
 2    continue
      ivnr=i1-1
      if (ivnr.lt.m) then
         print*,' !!!  Read Error within the Record: '
         print*,' You tried to read ',m,' variables '
         print*,' The following Record had only ',ivnr,' variables'
         goto 90
      endif

      return

 90   ierr=ios+1000
      print*,istr
      return 1

 99   ierr=ios
      return 1

      end





!gfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine gasflux(rad_lat,dday,dtair,dwind,dhumid,dsst, &
          dgfac,dca,direct,diffus,cosza,newLAI,xpsi,trans,pgros, &
          xpress,radsoil,cantemp,darkresp,cangl)
!c
!c     that's the 'Main'-Interface-Routine
!c--------------------------------------------------------------------------

      USE constants

      implicit none


      integer drow,dcol,dveg

      parameter (drow=6)
      parameter (dcol=18)
      parameter (dveg=1)

      common/res/transp,co2_fx,wdnlay,cantmp,cantotgl
      common/BigT/vegtyp
      common/rdsea/matrix(dcol,drow,dveg)
      common/pcontr/alat,date,live
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/vcontr/time,isirok
      common/ppsn2/o2,fvc,c,gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,fkc, &
          eavc,hdvc,dsvc,gmin,gmax,rdfac,ftau,eatau
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/vligt3/solan
      common/plight/lai,sai,angst

      double precision solan,cantotgl,xpsi
      double precision matrix,diffus,direct,cantemp,darkresp,cangl
      double precision dtair,dwind,dhumid,dsst,dgfac,dca
      double precision alat,date,live,transp,co2_fx,wdnlay, &
          cantmp,diff,dirn,tair,rh,sst,press,wind,ca,time
      double precision o2,fvc,c,gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,  &
          fkc,eavc,hdvc,dsvc,gmin,gmax,rdfac,ftau,eatau
      double precision ci,ga,pn,gs,gl,wd
      double precision par,tpairc,swtot,templf,therm,tran
      double precision lai,sai,angst,xpress
      double precision trans,pgros,newLAI,radsoil
      double precision cosza,rad_lat
      integer dday,vegtyp,veg
      integer isirok,j

      print*, 'function gasflux'


!cgy   only for common-block initializing
!cgy   solan in radians
      solan=Pi/2.-real(acos(cosza))
!
!c      if (solan .gt. 0.) then
!c        print*,solan,cosza,acos(cosza)
!c        read(*,*)
!c      endif

!c     latitude in radians
      alat=rad_lat

      diff=diffus
      dirn=direct

      vegtyp=1

      tair=dtair
      wind=dwind*100.    !m/s -> cm/s
      rh=dhumid/100.
      sst=dsst
      gfac=dgfac
      ca=dca

      tair=tair+273.16
      sst=sst+273.16

      date=dday

      lai=newLAI

      call initCommon

!cgy   changed, due to height dependent air pressure
!cgy   press reading is obsolete, should be changed
!cgy   xpress=press
      press=xpress

      call geom

      do j=1,2
         ci(j)=366.
         gs(j)=200.
      enddo

!c     calculate energy balance throughout canopy
      call baleng(radsoil)

!c     End-of-hour report
      call report(xpsi)

!c     Transpiration in l/m2h
      trans=transp
!c     CO2-Flux in umol/m2s
      pgros=co2_fx
!c     soilrespiration in mmol/m2h
!ccc      sresp=ssresp
      cantemp=cantmp
      darkresp=wdnlay
      cangl=cantotgl

      return
      end





!cgfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine initCommon

!c     Initialize the Common-Blocks for the site specific seasonal constants
!c--------------------------------------------------------------------------
      implicit none

      integer drow,dcol,dveg

      parameter (drow=6)
      parameter (dcol=18)
      parameter (dveg=1)

      common/rdsea/matrix(dcol,drow,dveg)
      common/BigT/vegtyp
      common/const/pi,rgas,stef
      common/ctrl1/needles,p_flag
      common/ctrl2/sigma,fac1
      common/pcontr/alat,date,live
      common/plight/lai,sai,angst
      common/fotpa/fop(22)
      common/penerg/ta,tr,tt,va,vr,vt
      common/ptemp/w
      common/pgeom/anglf
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/beps/omega

      double precision matrix,live,date,lai,alat,sai,omega, &
          angst,fop,ta,tr,tt,va,vr,vt,w, &
          anglf,pi,rgas,stef, &
          diff,dirn,tair,rh,sst,press,wind,ca,sigma
      double precision fac1,fac2,an,as,at
      integer i,a_flag,p_flag,vegtyp,needles


!c ----------------------------- control ----------------------------------
!c     Vegetationtype
      needles=int(matrix(1,1,vegtyp))
!c     Area basis of rates (see comments in subroutine psn6)
      a_flag=int(matrix(2,1,vegtyp))
!c     Leaf area basis (see comments in subroutine psn6)
      p_flag=int(matrix(3,1,vegtyp))
!e     maxLAI
!e     maxSAI
      omega=matrix(6,1,vegtyp)

!c     C,D are parameterized for Arbutus (Harley et al
!c     1986.  Oecologia 70:393-401)
!c     C: Scaling factor for co2-saturated photosynthesis
!c     D: Scaling factor for dark respiration (Ahrrenius function)
!c     ALPHA: Initial slope of the light response curve
!c     RDFAC: factor for DarkRespiration
!c     FOP(1....6,J1):---> FVC,C,D,O2,ALPHA,RDFAC

      do i=1,6
         fop(i)=matrix(i,2,vegtyp)
      enddo

!c     HA: Activation energy of Pml
!c     HD: Deactivation energy of Pml
!c     DS: Entropy factor
!c     EAKC: Activation energy for KC
!c     EAVC:  Activation energy for Vcmax
!c     FVC: Constant for Vcmax
!c     HDVC: Deactivation energy for Vcmax
!c     DSVC: Entropy term for Vcmax
!c     GMIN: Minimum conductance
!c     GMAX  MAXIMUM Conductance
!cEF   added EATAU, FTAU
!c     FOP(7...22,J1): ---> HA,HD,DS,E,EATAU,FTAU,EAKO,FKO,EAKC,FKC,
!c     EAVC,HDVC,DSVC,GMIN,GMAX,GFAC

      do i=7,22
         fop(i)=matrix(i-6,3,vegtyp)
      enddo

      anglf=matrix(1,4,vegtyp)
      live=matrix(2,4,vegtyp)
      w=matrix(3,4,vegtyp)
      angst=matrix(4,4,vegtyp)

!c     get fac1 and fac2 from matrix
      if (needles.eq.1) then
         fac1=matrix(1,5,vegtyp)
         fac2=matrix(2,5,vegtyp)
!c        Calculate sigma as in STANDFLUX
!c        Sigma is factor of bunching from Jarvis, James, Landsberg(1976),
!c        equation 18.
!c        Calculate sigma, where sigma=(An + At)/As, where An is total
!c        projected area of needles laid out separately, At is total
!c        projected area of twig without needles, As is total projected
!c        area of twig with needles.  2.57 converts total needle surface
!c        area to projected area for needles laid out separately.
!c        Calculations are done for case where projected area of intact
!c        twig is 1.0.
         an=1./(2.57*fac1)
         at=1/fac2
         as=1.
         sigma=(an+at)/as
      endif

!c     Calculate leaf area on a projected intact twig basis
!c     a_flag=0: lai is projected leaf area
!c     a_flag=1: lai is leaf surface area
      if (a_flag.eq.0) then
         if (needles.eq.1) lai=lai*2.57*fac1
      else
         if (needles.eq.1) then
            lai=lai*fac1
         else
!           for leaves "lai" is calculated for proj.leaf area
            lai=lai/2.
         endif
      endif


      return
      end






!gfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine geom

!     calculates static geometric attributes of canopy
!--------------------------------------------------------------------------
      implicit none

      common/const/pi,rgas,stef
      common/pcontr/alat,date,live
      common/plight/lai,sai,angst
      common/pgeom/anglf
      common/vgeom2/difac,difacu,difacb,difout,upvw,sivw,dnvw,ssfrac

      double precision fpass(2,9),favl(9),favls(9),vfac(9), &
          vfacu(9),zzi,zza,zzaa,zzii,zleaf,t1,t2,e1,e2,e3,aaa, &
          bbb,ccc,ddd,a2,b2,zz,zzz,ssup
      double precision pi,rgas,stef,alat,date,live,lai,sai, &
          angst,coslf,sinlf,anglf,penldf(9),pensdf(9),frac,  &
          verti,difac,difacu,difacb,difout,upvw,sivw,dnvw,ssfrac
      integer i

!     convert angles to radians; calculate commonly-used functions
!     of leaf angle, and convert leaf width to meters.

      angst=angst*pi/180.
      anglf=anglf*pi/180.
      coslf=cos(anglf)
      sinlf=sin(anglf)
!     w=w/100.

!     Calculate VFAC and VFACU arrays
!     Consider 9 annular bands, each representing a
!     portion of the sky as a source of diffuse radiation.
!     Ring 1 is on the horizon; ring 9 is at zenith.
      do 102 i=1,9
!        Calculate maximum (ZZA=90, 80, ... 10 degrees) and
!        minimum (ZZI = 80, 70, ... 0 degrees) zenith angles
!        that define this sky band.
         zza=(10-i)*10
         zzi=zza-10
!        Convert to radians
         zza=zza/180*pi
         zzi=zzi/180*pi
!c        The following calculations are based on Appendix B of
!c        Duncan et al. (1967).  They partition the leaf's view
!c        of the sky hemisphere into the fractions of the
!c        hemisphere that are potentially visible (assuming
!c        no obstacles from other leaves, etc) from
!c        upper (VFAC) and lower (VFACU) leaf surfaces.
!c        These view factors are computed numerically,
!c        one sky band at a time.

!c        The view factors are computed from the component
!c        terms E1, E2 and E3.  Each represents the fraction
!c        of the whole-sky hemisphere in a particular sky
!c        band.  E1 is calculated only for sky bands that
!c        are entirely visible from the leaf's top surface.
!c        E2 and E3 are calculated for sky bands visible
!c        from both the upper and lower surfaces.  E2 and E3
!c        represent the partitioning of a sky band to upper
!c        (E2) and lower (E3) surfaces.
!
!c        First, the sky band's angles are compared to the
!c        leaf's zenith angle to determine whether this band
!c        is visible only from the leaf's upper surface, or from
!c        both upper and lower (No sky band is visible only from
!c        the lower surface, since the sky band surrounds the
!c        leaf in all azimuthal directions).

!c        Zenith angle of leaf
         zleaf=pi/2.-anglf
         if (zzi.ge.zleaf) then
!c           This sky band is too low to be visible (in all
!c           azimuthal directions) from the top of the leaf.
!c           E1 won't be calculated.
            e1=0.
         else
!c           Calculate E1.  First, determine the angular limits.
            if (zzi.lt.zleaf .and. zza.gt.zleaf) then
!c              The leaf is within the sky band, so E1
!c              will be calculated for the portion of
!c              the sky band above the leaf.  Use a
!c              new lower sky band angle.
               zzaa=zleaf
            else
!c              The leaf is lower than the sky band. The entire
!c              sky band is visible only from the top surface.
               zzaa=zza
            endif
!c           Equation (11) of Appendix B in Duncan et al.: (1967)
            t1=sin(zzaa)
            t2=sin(zzi)
            e1=coslf*(t1*t1-t2*t2)
         endif
         if (zza.le.zleaf) then
!c           This sky band is visible in all azimuthal directions,
!c           and so it is inside the cone.  It doesn't contribute
!c           to E2 or E3.
            e2=0.0
            e3=0.0
         else
!c           Calculate E2 and E3.  First, determine anglular limits.
            if (zzi.lt.zleaf .and. zza.gt.zleaf) then
!c              The leaf is within the sky band, so E2
!c              and E3 will be calculated for the portion of
!c              the sky band below the leaf.  Use a
!c              new upper sky band angle.
               zzii=zleaf
            else
!c              The leaf is steeper than the sky band.
!c              Use the entire 10-degree sky band.
               zzii=zzi
            endif
!c           Partition view of sky in this band between top and
!c           bottom leaf surfaces:
!c           First, define sines and cosines whose angles
!c           will be found in calculation of E2 and E3 ...
            aaa=(1./tan(anglf))*(1./tan(zza))
            bbb=(1./tan(anglf))*(1./tan(zzii))
            ccc=(1./sinlf)*cos(zzii)
            ddd=(1./sinlf)*cos(zza)
!c           Limit AAA, BBB, CCC, DDD to 1.
            if (aaa.gt.1.) aaa=1.
            if (bbb.gt.1.) bbb=1.
            if (ccc.gt.1.) ccc=1.
            if (ddd.gt.1.) ddd=1.
            a2=-aaa
            b2=-bbb
!c           This band's contribution to the fraction of sky
!c           that is visible on top of leaf
!c           Equation (15) from Appendix B in Duncan et al (1967):
            e2=(coslf*((sin(zza)**2.)*acos(a2)-(sin(zzii)**2.) &
                *acos(b2))+asin(ccc)-asin(ddd)+cos(zzii)* &
                (abs(-cos(anglf+zzii)*cos(anglf-zzii)) &
                **(1./2.))-cos(zza)*(abs((-cos(anglf+zza) &
                *cos(anglf-zza)))**(1./2.)))/pi
!c           this band's contribution to the fraction of sky
!c           that is visible on bottom surface of leaf
!c           Equation (16) from Appendix B in Duncan et al (1967):
            e3=(coslf*((sin(zzii)**2.)*acos(bbb) &
                -(sin(zza)**2.)*acos(aaa))+asin(ccc)-asin(ddd) &
                +cos(zzii)*(abs((-cos(anglf+zzii)*cos(anglf &
                -zzii)))**(1./2.))-cos(zza)*(abs(-cos(anglf+zza) &
                *cos(anglf-zza))**(1./2.)))/pi
!c           Note: E1, E2 and E3 (all summed over the 9
!c           9 annular rings) should add up to 1.
         endif
!c        How much of the sky (in this sky band) is "visible"
!c        to side 1 of leaf in this layer?
         vfac(i)=e1+e2
!c        How much of the sky (in this sky band) is "visible"
!c        to side 2 of leaf in this layer?
         vfacu(i)=e3
 102  continue

!c     Consider 9 annular rings, each representing a
!c     portion of the sky as a source of diffuse radiation.
!c     Ring 1 is on the horizon; ring 9 is at zenith.
      do 104 i=1,9
!c        angle between horizon and middle of band, in
!c        degrees (5, 15, 25, ... 85)
         zz=(i-1)*10+5
!c        Convert to radians
         zz=zz/180*pi
!c        Calculate attenuation factors for diffuse radiation
!c        due to leaves (PENLDF) and stems (PENSDF).
!c        These coefficients are analaogus to ATTLEV and ATTSTM
!c        in the direct light calculations.  They are ratios
!c        of the leaf's (or stem's) shadow to its area.  The
!c        ratio is 1 for a horizontal leaf and overhead light
!c        source; it is 0 when the light source angle equals
!c        the leaf angle.  Reference: Appendix A of Duncan et al.
!c        (1965), based on the work of Wilson and Reeve.
         call penet(anglf,zz,lai,zzz,penldf(i))
         call penet(angst,zz,sai,zzz,pensdf(i))
 104  continue

!c     FRAC and VERTI only also depend only on leaf geometry.
!c     Calculate an interception factor for vertically-moving radiation.
!c     This is simply a measure of projected leaf area.
      verti=lai*coslf
      if (verti.gt.1.) then
         verti=1.0
      endif

!c     Calculate fractional diffuse radiation from leaf within a
!c     "vertical cone," using Duncan et al.'s equation (2)
!c     where Beta =zenith angle .57 radians.  The fraction is
!c     1 for horizontal leaves, and decreases as angle increases.
!c     frac=pi*coslf*(.573576)**2.
!c     if (frac.lt.0.5) then
!c        frac=(0.5-coslf)/2.0+coslf
!c     elseif (frac.gt.1.0) then
!c        frac=1.0
!c     endif
!c     set frac to (1-cos(35 degrees))
      frac=0.18085

      upvw=0.0
      sivw=0.0
      difac=0.0
      difacu=0.0
      difacb=0.0

      difout=0.
!c     Calculate fraction of the available diffuse ir or sw in the
!c     sky band that will pass unimpeded into the layer.
      do 108 i=1,9
!c        All of the above-canopy radiation enters top layer.
         fpass(1,i)=1.
!c        Fraction passing into next layer down depends on
!c        penetration through leaves and stems in this layer.
         fpass(2,i)=fpass(1,i)*penldf(i)*pensdf(i)
!c        Fraction of sky diffuse light avail. to leaves and stems
!c        is the average of that entering and exiting the layer.
         favls(i)=(fpass(1,i)+fpass(2,i))/2.
!c        Fraction available to leaves only.
         favl(i)=favls(i)*(1-penldf(i))/   &
             ((1-penldf(i))+1-pensdf(i))

!c        Calculate maximum (ZAA=90, 80, ...10 degrees) and
!c        minimum (ZII = 80, 70,...0 degrees) zenith angles
!c        that define the sky band.
         zza=(10-i)*10
         zzi=zza-10
!c        convert to radius
         zza=zza/180*pi
         zzi=zzi/180*pi
!c        Accumulate fraction of of sky diffuse incident on soil.
         difout=difout+(sin(zza)**2-sin(zzi)**2)*fpass(2,i)

!c        Calculate diffuse factors for top leaf surface (DIFAC),
!c        under surface (DIFACU), and both surfaces (DIFACB).
!c        Each factor is calculated as the sum over sky bands of
!c        the products of the diffuse fractions and the view factors.
!c        The factors represent the fraction of the above-canopy
!c        diffuse irradiance intercepted by leaves.
         difac=difac+favl(i)*vfac(i)
         difacu=difacu+favl(i)*vfacu(i)
         difacb=difacb+favl(i)*(vfac(i)+vfacu(i))
         if (i.gt.4) then
!c           Calculate composite longwave view factor ("up view")
!c           for interception of leaf diffuse ir from the
!c           layer above.
            upvw=upvw+(vfac(i)+vfacu(i))*(1-fpass(1,i))
         else
!c           Calculate composite longwave view factor ("side view")
!c           for interception  of leaf diffuse ir emitted from the
!c           same layer.
            sivw=sivw+(vfac(i)+vfacu(i))*(1-fpass(1,i))
         endif
 108  continue

!c     All soil-surface emitted ir makes it up into bottom layer.
      ssup=1.
!c     Fraction of soil-surface emitted ir intercepted in layer
      ssfrac=ssup*coslf
!c     Calculate special longwave view factor ("down view")
!c     for interception  of leaf diffuse ir emitted from the
!c     layer below.
      dnvw=1-ssfrac
!c     Fraction of soil-surface-emitted IR getting into next layer.
      ssup=ssup*exp(-verti)

      return
      end





!cgfx+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine penet(incl,angle,oai,ff,g)
!c
!c     calculates radiation penetration through stems or leaves
!c
!c     argument list:
!c       incl   -- mean inclination angle for stem or leaf
!c       angle  -- angle from horizon of source of light
!c       oai    -- area index (leaf or stem area index)
!c       ff     -- the Wilson-Reeve ratio calculated here.
!c       g      -- the fraction of light passing thru a layer
!c                 composed only of leaves or stems.
!c--------------------------------------------------------------------------
      implicit none

      double precision incl,angle,oai,ff,g,sinlit,coslit,sininc, &
          cosinc,a1,th,p

      sinlit=sin(angle)
      coslit=cos(angle)
      sininc=sin(incl)
      cosinc=cos(incl)

      if (angle.gt.0.001) then
!c        Source of light is above horizon.  Proceed with calculations.
!c        Calculate FF, the Wilson-Reeve ratio F'/F where
!c        F' is the area of the shadow cast by a leaf or stem
!c        on a plane normal to the irradiation;
!c        F  is the actual area of the leaf.
!c        As derived in Duncan et al. 1967. Hilgardia 38:181-205.
         if (incl.gt.angle) then
!c           the lower surface of the leaf or stem is partially
!c           illuminated.  Complex case.
            a1=acos(cosinc/sininc*sinlit/coslit)
            th=(1.-a1/1.570796)*cosinc*sinlit
            ff=.6366198*sininc*coslit*sin(a1)+th
         else
!c           The lower surface cannot be illuminated.  Simple case.
            ff=cosinc*sinlit
         endif

!c        Calcuate equation (3) of Duncan et al (1967) for P,
!c        which is supposedly the probability of a photon touching
!c        a single leaf (actually P*exp(-P) is that probablity)
         p=oai*ff/sinlit

!c        Calculate the fraction of light passing through
!c        a canopy layer (composed only of leaves or stems)
!c        Source: Duncan et al. (1967) p.188
         g=exp(-p)
      else
!c        Source of light is at or below horizon.  Don't bother.
         g=0
      endif

      return
      end





!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine baleng(radsoil)
!c
!c     calculates energy balance throughout canopy
!c-----------------------------------------------------------------------
!c                         Output Variables Dictionary
!c-----------------------------------------------------------------------
!c       Name                   Description                       Units
!c
!c     EWAIR        Vapor pressure of air above canopy             mB
!c     ISIROK       Flag indicating whether sky and soil          none
!c                    longwave have been solved
!c     J            Index for sun(1) or shade(2)                  none
!c     L            Index for layer number                        none
!c     TEMPLF(L,J)  Leaf temperature                              deg C
!c     TLFK         Leaf temperature                              deg K
!c     TPAIR(L)     Temperature of air in layer L                 deg K
!c     TPAIRC(l)    Temperature of air in layer L                 deg C
!c     WIND(L)      Wind velocity above canopy(later attenuated   cm s-1
!c                   through canopy)
!c-----------------------------------------------------------------------
      implicit none

      common/pcontr/alat,date,live
      common/vcontr/time,isirok
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/vtemp/tpair,ewair,ewleaf(2)
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)
      common/vligt1/aldif,ttsi,sarea
      common/vligt3/solan
      common/plight/lai,sai,angst

      double precision aldif,ttsi,sarea
      double precision lai,sai,angst
      double precision eloss(2),esat,tlfk
      double precision alat,date,live,time,diff,dirn,tair,rh, &
          sst,press,wind,ca,tpair,ewair,ewleaf, &
          par,tpairc,swtot,templf,therm,tran
      double precision ci,ga,pn,gs,gl,solan,wd,humair,radsoil
      integer isirok,j

!c     Canopy Light Interception
      call shrtwv(lai,sai,radsoil)

      if (wind.le.3.) wind=3.

!c     Saturated vapor pressure of air above canopy
      call vap(tair-273.16,ewair)

!c     Actual vapor pressure
      ewair=ewair*rh

!c     Air Temperature Profile
!c     Calculate a logarithmic air temperature profile in the canopy from
!c     air temperature above the canopy and the soil surface temperature
!c     if there are NO PROFILE VALUES
      tpair=tair

!c     Calculate temperature of air in degrees C (for use in REPORT)
      tpairc=tpair-273.16
!c     The first approximation of leaf temperatures for both sunlit
!c     and shaded leaves is the air temperature.
      templf(1)=tpairc
      templf(2)=tpairc

!c     Calculate Relative Humidity of air in layer if there are
!c     NO PROF. VALUES
!c     First determine saturated v.p. at this temperature.
      call vap(tpairc,esat)
!c     R.H. is the ratio of the v.p.s
      humair=ewair/esat         ! never used

!c     Canopy IR Budget and Leaf Temperature Iterations
!c     Set the control variable ISIROK to 0 so the hour's first call
!c     to LONGWV will include calculations of DIFSIR(), UPVW() and SIVW().
      isirok=0

!c     Pass through the canopy from top to bottom, for a second
!c     approximation of leaf temperatures as affected by longwave
!c     radiation from adjacent leaves and by transpiration.
!c     Calculate leaf interception of longwave radiation
      call longwv
!c     Make a second approximation of temperature of sunlit leaves.
      call ltemp(1)
!c     Make a second approximation of temperature of shaded leaves.
      call ltemp(2)

!c     Pass through the canopy, from bottom to top, for a third
!c     approximation of leaf temperatures.
!c     Recalculate leaf interception of longwave, taking into account
!c     the longwave re-radiation at the second-approximation
!c     leaf temperatures.
      call longwv
!c     Make a third approximation of temperature of sunlit leaves.
      call ltemp(1)
!c     Make a third approximation of temperature of shaded leaves.
      call ltemp(2)

!c     Canopy Conductance and Transpiration
!c     Pass through the canopy from the top down, to calculate
!c     stomatal conductance and transpiration as a function of
!c     the equilibrium leaf temperatures.
!c     Consider sunlit (J=1 and shaded (J=2) leaves
      do 26 j=1,2
!c        Leaf temperature, absolute
         tlfk=templf(j)+273.16
!c        Calculate stomatal conductance and transpiration;
!c        keep track of energy loss
         call htloss(tlfk,eloss(j),j,1)
 26   continue

      return
      end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine shrtwv(lai,sai,radsoil)
!c
!c     Calculate canopy light interception.
!c
!c     This subroutine is called by the subroutine baleng
!c-----------------------------------------------------------------------
!c                         Output Variables Dictionary
!c-----------------------------------------------------------------------
!c       Name                   Description                       Units
!c
!c     PAR(L,J)     Photosynthetic photon flux density          u mol m-2
!c     R(L)         Reflectance for PAR or total shortwave      proportion
!c     SWTOT(L,J)   Total shortwave radiation reaching leaves    kW m-2
!c     T(L)         Transmittance for PAR or total shortwave    proportion
!c-----------------------------------------------------------------------
      USE constants
      implicit none

      common/pcontr/alat,date,live
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/vligt1/aldif,ttsi,sarea
      common/vligt2/difs,difs2
      common/vligt3/solan
      common/vgeom2/difac,difacu,difacb,difout,upvw,sivw,dnvw,ssfrac
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/beps/omega


      double precision alat,date,live,diff,dirn,tair,rh,sst,press, &
          wind,ca,aldif,ttsi,sarea,difs,difs2,solan, &
          difac,difacu,difacb,difout,upvw, &
          sivw,dnvw,ssfrac,par,tpairc,swtot,templf,therm,tran
      double precision dir_in,theta_in,omega,lai,sai,theta_avg,  &
          sdif_under,sshade,ssun,radsoil


!c     calculate solar altitude and irradiance from sun and sky
!cgy   solan,diff,dirn is calculated above
!cgy      call sunsky

!c     calculate direct radiation on horizontal surface
!ce     dir_in=dirn*sin(solan)   ! normal radiation -> horizontal
      dir_in=dirn               ! dirn is already horizontal

!c     calculate solar zenith angle
      theta_in=pi/2.-solan

      sarea=0.
      aldif=0.
      ttsi=0.
      swtot(1)=0.
      swtot(2)=0.
      par(1)=0.
      par(2)=0.
      if (solan.gt.0.) then
!c        average sunlit leaves
!c        SAREA Area index of leaves irradiated directly by sun
         sarea=2.*cos(theta_in)*(1.-exp(-0.5*omega*lai/   &
             cos(theta_in)))
!c        Correct for effect of stems. I have not used lai+sai in the
!c        above line, but did the following correction. Might be worth
!c        trying lai+sai above and remove the following line.
         sarea=sarea*lai/(lai+sai)

!c        radiation from multiple scattering (assume that SAI does not scatter,
!c        and used only LAI here!)
!c        ALDIF     Total leaf diffuse shortwave intercepted by leaves
         aldif=0.07*omega*dir_in*(1.1-0.1*lai)*exp(-cos(theta_in))

!c        sdif_under
         theta_avg=0.537+0.025*(lai+sai)
         sdif_under=diff*exp(-0.5*omega*(lai+sai)/theta_avg)

!c        diffuse sky radiation
         sshade=(diff-sdif_under)/(lai+sai)
!c        direct sun radiation
         ssun=0.5*dir_in/cos(theta_in)

!c        TTSI Shortwave irradiation from sun and sky diffuse reaching leaves
         ttsi=ssun+sshade

!c        SWTOT(J)   Total shortwave radiation reaching leaves
!c        Calculate total sw in sunlit leaves for use in LTEMP
         swtot(1)=aldif+ttsi

!c        Calculate sw in shade
         swtot(2)=aldif+sshade

!c        PAR(J)     Photosynthetic photon flux density
!c        Convert kilowatts/m2 to micromoles/m2/sec using empirical factor
!c        !!! This however assumes that reflectance and transmittance are
!c        the same for PAR and total shortwave - in contrast to original
!c        GASFLUX !!!
!c        Sunlit leaves
         par(1)=swtot(1)*2.24e03
!c        Shaded leaves
         par(2)=swtot(2)*2.24e03
!c        What reaches the soil?
         radsoil=dir_in+diff-swtot(1)*sarea-swtot(2)*(lai+sai-sarea)
      endif

      return
      end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine vap(tc,ew)
!c
!c     calculate saturated water vapor pressure
!c
!c     this subroutine is called every hour by the main program, and
!c     several times by baleng.
!c
!c     argument list:
!c      temp  -- temperature of air or leaf (degrees C)
!c      ew    -- partial pressure of water vapor (saturated v.p., in mB).
!c-----------------------------------------------------------------------
      implicit none

      double precision tc,ew

!c     Use either the mechanistic or the empirical function
!c     Mechanistic equation
!c     Source: List (1984) Smithsonian Meteorological Tables, p.350.
!c     tk=tc+273.16
!c     x1=11.344*(1.-tk/373.16)
!c     xxx=373.16/tk
!c     xx1=-3.49149*(xxx-1.)
!c     t1=exp(x1*log(10.))
!c     t2=exp(xx1*log(10.))
!c     xogen=-7.90298*(xxx-1.)+5.02808*alog10(xxx)-
!c     +      1.3816*1.e-07*(t1-1.)+8.1328*.001*(t2-1.)+3.0057149
!c     The next line is the equivalent of  EW=10.**XOGEN
!c     (2.302 is natural log of 10).
!c     ew=exp(xogen*2.30258509)
!c
!c     To speed up VAP, comment-out the above lines and un-comment the one
!c     below that calculates EW.  This new function is from Teten, and
!c     is described in  Beers, N.R.  1945.  Meteorolgical thermodynamics
!c     and atmospheric statics.  IN: F.A. Berry, E. Bollay and N.R. Beers
!c     (eds), Handbook of Meteorology.  Section V.  pp. 313-409.  McGraw-Hill.

      ew=6.11*exp(17.26938818*tc/(237.3+tc))

      return
      end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine longwv
!c
!c     calculates leaf interception of longwave (ir) radiation from
!c     sky, leaves and soil
!c
!c     This subroutine is called by baleng.
!c-----------------------------------------------------------------------
!c                         Output Variables Dictionary
!c-----------------------------------------------------------------------
!c       Name                   Description                       Units
!c
!c     DOWNIR(L)    Total longwave emitted by leaves in layer     kW m-2
!c                    L-1 intercepted in layer L
!c     THERM(L)     Total longwave radiation reaching leaves      kW m-2
!c     UPIR(L)      Total longwave emitted by leaves in layer
!c                    L+1 intercepted in layer L                  kW m-2
!c-----------------------------------------------------------------------
      implicit none

      common/const/pi,rgas,stef
      common/pcontr/alat,date,live
      common/vcontr/time,isirok
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/plight/lai,sai,angst
      common/vligt1/aldif,ttsi,sarea
      common/vgeom2/difac,difacu,difacb,difout,upvw,sivw,dnvw,ssfrac
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/vtemp/tpair,ewair,ewleaf(2)

      double precision ir,lai,live,pi,rgas,stef,alat,date,time,diff,   &
            dirn,tair,rh,sst,press,wind,ca,sai,angst,  &
            aldif,ttsi,sarea,difac,difacu,difacb,difout,  &
            upvw,sivw,dnvw,ssfrac,par,tpairc,swtot,templf,  &
            therm,tran,tpair,ewair,ewleaf,ea,sst4,bb1,bb2,t1
      double precision resir,eirlv,air,difsir,soilir
      integer isirok

      save difsir,soilir


!c     The following calculations need to be done once per hour.
!c     They are skipped if the flag ISIROK has been reset to 1.

      if (isirok.eq.0) then
!c        Longwave Radiation from Sky
!c        Above-canopy longwave (Brutsaert 1975)
         ea=.58*((.77586*ewair)**(1./7.))
         ir=ea*stef*(tair**4.)
!c        Track the sky longwave through the canopy.
!c        DIFSIR is the total sky diffuse ir received by leaves.
!c        DIFACB, calculated in GEOM, is a function of the sky ir,
!c        the attenuation of sky ir through the canopy, and
!c        "view factors" for top and bottom leaf surfaces.
          difsir=ir*difacb

!c        Longwave Radiation from Soil
!c        How much IR is the soil radiating?
         sst4=(sst**4)*stef
!c        How much of that IR is intercepted by leaves in each layer?
!c        Interception of soil IR by leaves (per leaf area)
         soilir=sst4*ssfrac
!c        Reset the flag ISIROK.
         isirok=1
      endif

!c     Leaf-emitted Longwave Radiation
!c     Calculate emitted infrared as a weighted average of
!c     the emitted IR from sunlit and shaded leaves.
      bb1=templf(1)+273.16
      bb1=bb1*bb1
      bb1=bb1*bb1
      bb2=templf(2)+273.16
      bb2=bb2*bb2
      bb2=bb2*bb2
!clai
      if (lai.gt.0.) then
         t1=sarea/lai
         eirlv=(t1*bb1+(1.-t1)*bb2)*stef
      else
         eirlv=0.
      endif
!c     Residual IR is emitted from leaves but stays in the layer.
      resir=sivw*eirlv

!c     Summarize IR interception.
!c     Sum of leaf-emitted IR reaching leaves in this layer.
      air=resir
!c     Total ir reaching leaves in layer
      therm=difsir+air+soilir

      return
      end





!
!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine ltemp(j)
!c
!c     iteratively calculates equilibrium leaf temperature
!c
!c     This subroutine is called by the baleng subroutine.
!c
!c     argument list:
!c       l   -- number of canopy layer
!c       j   -- sunlit (j=1) or shaded (j=2) leaves.
!c-----------------------------------------------------------------------
      implicit none

      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)

      double precision ai,ta,par,tpairc,swtot,templf,  &
          therm,tran,absrad,tl1,tl2,tl3,diss1,diss2,diss3,   &
          slope,zdif,errval,dif
      integer j,maxit,iter

!cgy   absorptance of total shortwave
      ta=0.5
!cgy   absorptance of infrared
      ai=0.975

!c     Determine the equilibrium leaf temperature, i.e., the
!c     temperature at which energy losses equal energy gains.
!c     Total absorbed energy.  A leaf temperature will be found
!c     where energy losses equal this value.
      absrad=ta*swtot(j)+ai*therm

!c     When air temperatures are too low, because there is no convective
!c     heat losses, leaf temperatures can be too far below air temperatures
!c     In thoses cases we set leaf temperature equal to air temp at aech layer
!c      if (tpairc.le.10.) then
!c         templf(j)=tpairc
!c         goto 2
!c      endif

!c     Calculate the heat loss at the current estimate of leaf
!c     temperature.
      tl1=templf(j)+273.16

!c     Radiation gains and losses will differ by no more than ERRVAL
      errval=0.005

      call htloss(tl1,diss1,j,1)

!c     Is the initial estimate good enough?
      dif=diss1-absrad
      if (abs(dif).ge.errval) then
!c        An iterative solution will be necessary.  Solve for
!c        equilibrium by means of iterative linear interpolation:
!c        make an overestimate and an underestimate of leaf temperature,
!c        then calculate energy losses for each estimate.  Calculate an
!c        intermediate temperature by assuming linearity between the two
!c        points and interpolating for the temperature with the given
!c        heat loss.  Determine the difference between this and the
!c        equilibrium heat loss.  If this difference is large
!c        enough, set either the upper or lower temperature bound
!c        to the intermediate one, and repeat.  Since the relation
!c        between heat loss and leaf temperature is nearly linear,
!c        two iterations are almost always sufficient.

!c        Maximum number of iterations
         maxit=10
         iter=0

!c        Set upper and lower temperature bounds for iteration
!c        The initial estimate's use as upper or lower limit
!c        depends on whether the initial estimate is greater or
!c        smaller than the total absorbed energy.
         if (diss1.gt.absrad) then
!c           Use initial estimate as upper bound
            tl2=tl1
            tl3=templf(j)+263.16
            diss2=diss1
            call htloss(tl3,diss3,j,2)
         else
!c           Use initial estimate as lower bound
            tl2=templf(j)+283.16
            tl3=tl1
            diss3=diss1
            call htloss(tl2,diss2,j,2)
         endif

!c        Start the iteration loop
 1       continue
         iter=iter+1
!c        Calculate an intermediate temperature, using the point-slope
!c        form of the linear equation TL=a+b*DISS.
         zdif=diss2-diss3
         slope=(tl2-tl3)/zdif
         tl1=tl2+slope*(absrad-diss2)
!c        Heat loss at this intermediate temperature.
         call htloss(tl1,diss1,j,3)
!c        Close enough?
         dif=diss1-absrad
         if ((iter.le.maxit).and.(abs(dif).ge.errval)) then
!c           Set either the upper or lower bound to the inter.temperature
            if (dif.gt.0.0) then
!c              equilibrium temp. must be below tl3; reduce upper bound
               tl2=tl1
               diss2=diss1
            elseif (dif.lt.0.0) then
!c              equilibrium temp. must be above tl3;
!c              increase lower bound
               tl3=tl1
               diss3=diss1
            endif
!c           Continue the iteration
            goto 1
         endif
      endif

!c     Set leaf temperature to the final mid-temperature.
      templf(j)=tl1-273.16

 2    return
      end



!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine htloss(tk,d,j,iflag)
!
!c     calculates leaf gross energy loss at a given temperature
!c
!c     This subroutine is called by the baleng or the ltemp subroutine.
!c
!c     argument list:
!c       tk  -- temperature of leaf, Kelvin
!c       d   -- energy loss from leaf at this temperature
!c       l   -- number of canopy layer
!c       j   -- sunlit (j=1) or shaded (j=2) leaves.
!c       iflag -- time of call.  1=first time for leaf; 2=second.
!c-----------------------------------------------------------------------
!c     Data Dictionary
!c     ===============
!c     1. Subroutine Inputs
!c     Name         Units       Description
!c     AI(L)        proportion  Leaf absorptance of longwave
!c     EWAIR        mB          Vapor pressure of air
!c     EWLEAF(L,J)  mB          Vapor pressure of leaf air space
!c     PRESS        mB          Atmospheric pressure
!c     RGAS                     The gas constant
!c     STEF                     Stefan-Bolzmann constant
!c     TK           degrees K   Current estimate of this leaf's temperature
!c     TPAIR(L)     degrees K   Temperature of air in layer
!c     W(L)         cm          Leaf width, or characteristic dimension
!c     WIND(L)      cm s-1      Wind speed in layer
!c
!c     2. Subroutine Outputs
!c     Name         Units       Description
!c     CON(L,J)     kW m-2      Convective heat loss
!c     D            kW m-2      Energy lost from leaf
!c     ELAT(L,J)    kW m-2      Latent heat of evaporation
!c     EMITL(L,J)   kW m-2      Heat dissipation due to longwave
!c                              radiant emittance
!c     GA           mmol m-2 s-1 Boundary layer conductance of leaf
!c     HUMBL(L,J)   proportion  Relative humidity within the boundary
!c                              layer
!c     TRAN(L,J)    mmol m-2 s-1 Transpiration from leaf
!c
!c     3. Internal Variables
!c     Name         Units       Description
!c     DBL          m           Thickness of boundary layer
!c     DWV          m2 s-1      Diffusion coefficient for water
!c     RTK                      Gas constant * Leaf temperature
!c     SQWID                    Square root of leaf width
!c     SQWIN                    Square root of wind speed
!c-----------------------------------------------------------------------

      USE constants
      implicit none

      common/ctrl1/needles,p_flag
      common/ctrl2/sigma,fac1
      common/pcontr/alat,date,live
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/ptemp/w
      common/vtemp/tpair,ewair,ewleaf(2)
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)
      common/vligt1/aldif,ttsi,sarea

      save dbl
      double precision live,kair,alat,date,diff,dirn, &
            tair,rh,sst,press,wind,ca,w,tpair,ewair,  &
            humbl(2),ewleaf,elat(2),emitl(2),par,fac1, &
            swtot,templf,therm,tran,ci,ga,pn,gs,gl,tpairc, &
            tk,d,tcbl,dbl,sqwid,sqwin,tc,rtk,dwv,elatnt,sigma, &
            con(2),ai
      double precision aldif,ttsi,sarea,wd
      integer j,iflag,needles,p_flag

      if (iflag.eq.1) then
!c        First time through this subroutine for the current leaf.
!c        Square roots of leaf width and wind
         sqwid=w**.5
         sqwin=wind**.5
!c        Thickness of boundary layer (in m).  Nobel 1983, eq. 7.7
!c        W and Wind should both be meters or centimeters
         dbl=0.004*sqwid/sqwin
      endif

      tc=tk-273.16
      rtk=rgas*tk
!c     Temperature of boundary layer, used in temp. sensitivities
      tcbl=(tc+tpairc)/2.

!c     The three components of the leaf's heat losses to the environment are
!c     1. Radiation losses,
!c     2. Conductive/convective losses -- also known as sensible heat,
!c     3. Latent heat of evaporation.
!c     All of these heat losses are temperature dependent.
!c     Consider each component in turn.

!c     ------------------------------------------------------------------
!c     1. Radiation losses
!c     Heat dissipation due to longwave radiant emittance
!c     Assumes that emissivity and absorptance coefficients are equal.
!cgy   absorptance of infrared
      ai=0.975
      emitl(j)=tk*tk*tk*tk*stef*ai*2.
!c     ------------------------------------------------------------------

!c     2. Conductive/convective losses
!c     Thermal conductivity coefficient of air.
!c     Temperature sensitivity calibrated from Nobel (1983) Appendix II
      kair=.024343+6.7143e-05*tcbl
!c     Convective heat loss: Nobel (1983) eq. 7.10 (kW m-2)
      con(j)=2.*kair/dbl*(tk-tpair)/1000.
!c     ------------------------------------------------------------------

!c     3. Latent heat of evaporation.  Requires calculation of transpiration.
!c     (and photosynthesis if it is calculated in conjunction with
!c     stomatal conductance)

!c     Diffusion coefficient for water vapor (m2 s-1)
!c     Temperature sensitivity calibrated from Nobel (1983) Appendix II
!c     dwv=2.126e-05 + 9.2e-08 * tc  / old statement prob. adapt. for CO2
      dwv=2.126e-05+1.48e-07*tc
!c     Boundary layer conductance to H2O: Nobel (1983) eq. 8.3 (mm s-1)
      ga=dwv/dbl*1000.
!c     Boundary layer conductance to H2O: from Nobel (1983) eq. 8.8.
!c     mmol m-2 s-1
      ga=ga*(press*100.)/rtk

!c     Modify boundary layer for bunching of needles:
!c     from Jarvis,James,Landsberg(1976) eq. 20
      if (needles.eq.1) ga=ga/(1.67*sigma**0.43)

!c     Calculate water vapor pressure of leaf based on
!c     this temperature, assuming humidity is saturated
      call vap(tc,ewleaf(j))
!c     Estimate boundary layer humidity.
      humbl(j)=ewair/ewleaf(j)

      if (humbl(j).gt.1.) humbl(j)=1.
!c     Stomatal conductance (GS).
!c     CHANGE: Insertion of the subroutine ASSPAR which assigns the model
!c     parameters for each species according to each layer and each species
      call asspar

      call psn6(j,par(j),tk,humbl(j))

!c     Total leaf conductance to H2O: mmol m-2 s-1
      if (gs(j).le.0.) gs(j)=0.00001
      gl(j)=1./(1./ga+1./gs(j))

!c     Transpiration in mmol m-2 s-1
      tran(j)=(ewleaf(j)-ewair)/press*gl(j)

!c     Prevent transpiration from going negative.
      if (tran(j).lt.0.) tran(j)=0.

!c     Temperature function for latent heat
!c     of vaporization (Kjoules/mmol water)
      elatnt=.0450643-4.29143e-05*tc
!c     latent heat loss kW/m2
      elat(j)=elatnt*tran(j)
!c     ------------------------------------------------------------------
!c     Total energy loss (kilowatts per square meter)
      d=emitl(j)+elat(j)+con(j)

      return
      end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine asspar
!c
!c     This subroutine assigns the model parameter values from the commonblock
!c     /FOTPA/FOP(30,ISPEC) to the values which are needed in the simulation
!c     subroutine PSN6
!c     ISPEC: Array with the species of the various layers  (NLAY)
!c     ISP:   Index of the species
!c     L:     Index of the canopy layer
!c-----------------------------------------------------------------------
      implicit none

      common/BigT/vegtyp
      common/FVCcalc/seasonstart,seasonend
      common/tttest/cutcntt
      common/fotpa/fop(22)
      common/ppsn2/o2,fvc,c,gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,fkc,   &
          eavc,hdvc,dsvc,gmin,gmax,rdfac,ftau,eatau
      common/pcontr/alat,date,live

      double precision live,fop,alat,date,o2,fvc,c,rdfac,ftau,eatau,  &
          gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,fkc,eavc,hdvc,dsvc, &
          gmin,gmax,FVCstart,FVCend
      integer cutcntt,vegtyp,seasonstart,seasonend

!c     C,D are parameterized for Arbutus (Harley et al
!c     1986.  Oecologia 70:393-401)
!c     C: Scaling factor for co2-saturated photosynthesis
!c     D: Scaling factor for dark respiration (Ahrrenius function)
!c     ALPHA: Initial slope of the light response curve
!c     RDFAC: factor for dark respiration
!c     FOP(1....6,n):---> FVC,C,D,O2,ALPHA,RDFAC
!
!c     HA: Activation energy of Pml
!c     HD: Deactivation energy of Pml
!c     DS: Entropy factor
!c     EAKC: Activation energy for KC
!c     EAVC:  Activation energy for Vcmax
!c     FVC: Constant for Vcmax
!c     HDVC: Deactivation energy for Vcmax
!c     DSVC: Entropy term for Vcmax
!c     GMIN: Minimum conductance
!c     FOP(7...20,n): ---> HA,HD,DS,E,EATAU,FTAU,EAKO,FKO,EAKC,FKC,EAVC,HDVC,
!c     DSVC,GMIN,GMAX

      fvc=fop(1)
      c=fop(2)
      d=fop(3)
      o2=fop(4)
      alpha=fop(5)
      rdfac=fop(6)
      ha=fop(7)
      hd=fop(8)
      ds=fop(9)
      e=fop(10)
      eatau=fop(11)
      ftau=fop(12)
      eako=fop(13)
      fko=fop(14)
      eakc=fop(15)
      fkc=fop(16)
      eavc=fop(17)
      hdvc=fop(18)
      dsvc=fop(19)
      gmin=fop(20)
      gmax=fop(21)
!ce    gfac=fop(22)

      if (vegtyp.eq.1) then               !rice fvc=70.
         fvc=-0.5*date+170.                   !jd200=70 jd280=30
         if (fvc.ge.70.) fvc=70.
         if (fvc.lt.0.) fvc=0.
      endif

      alpha=0.0008*fvc
      if (alpha.gt..06) alpha=0.06
      c=fvc/2.1
      d=fvc*0.025

      return
      end







!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine psn6(j,phar,tk,humid)
!c
!c     Leaf Photosynthesis Simulation
!c     Farquhar/Wong model of ps and g.
!c     Stomatal conductance: Ball et al. 1987
!c     ECOCRAFT/LTEEF equations
!c
!c     p_flag: Flaechenbasis der Raten, die durch PSParam berechnet werden
!c             0 falls Raten pro LAI, 1 falls Raten pro Oberflaeche
!c             berechnet werden
!c     a_flag: (fuer Bestandesmodell)
!c             Art der Blattflaeche: 0 falls Blattflaeche in LAI
!c             1 falls Blattflaeche als Leaf surface area angegeben wird
!c-----------------------------------------------------------------------
      USE constants
      implicit none

      common/ctrl1/needles,p_flag
      common/ctrl2/sigma,fac1
      common/pcontr/alat,date,live
      common/vdriv/diff,dirn,tair,rh,sst,press,wind,ca
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)
      common/ppsn2/o2,fvc,c,gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,fkc,  &
          eavc,hdvc,dsvc,gmin,gmax,rdfac,ftau,eatau
      common/vtemp/tpair,ewair,ewleaf(2)

      double precision live,alat,date,diff,dirn,  &
          tair,rh,sst,press,wind,ca,ci,ga,pn,wd,gs,gl, &
          o2,fvc,c,gfac,d,alpha,ha,hd,ds,e,eako,fko,eakc,fkc,  &
          eavc,hdvc,dsvc,gmin,gmax,tpair,ewair,fac1, &
          ewleaf,phar,tk,rtk,vcm,dd,pml,pm,rd, &
          ftau,eatau,tau,gamma,cint,pnet,sigma,cint1

      double precision kc,ko,fac,aaa,bbb,ddd,eee,wc,wj,gpres,rdfac,humid
      integer j,needles,flag,p_flag


!c     Gas constant * Absolute temperature
      rtk=rgas*tk

!c     Farquar and Wong (1984, Aust. J. Plant Phys. 11:191-210)
!c     Michaelis constant for O2
      ko=fko*exp((eako*(tk-298.))/(rtk*298.))
!c     Michaelis constant for CO2
      kc=fkc*exp((eakc*(tk-298.))/(rtk*298.))

!c     Maximum rate of RuP2 carboxylation
      vcm=fvc*exp((eavc*(tk-298.))/(rtk*298.))/  &
          (1+exp((dsvc*tk-hdvc)/rtk))* &
          (1+exp((dsvc*298.-hdvc)/(8.31*298.)))

!c     Calculate PML, the assimilation rate at saturating
!c     PAR and saturating co2, i.e. Jmax/4.
!c     Johnson et al. (1942. J Cell Comp Physiol 20:247-268) equation
      pml=c*exp((ha*(tk-298.))/(rtk*298.))/   &
          (1+exp((ds*tk-hd)/rtk))*  &
          (1+exp((ds*298.-hd)/(8.31*298.)))

!c     Smith equation
      dd=1+alpha*alpha*phar*phar/(pml*pml)
      pm=alpha*phar/sqrt(dd)

      rd=d*exp((e*(tk-298.))/(rtk*298.))

      tau=ftau*exp((eatau*(tk-298.))/(rtk*298.))
      gamma=.5*o2*1000./tau

!c     if (phar.le.0.) fac=1
      if (phar.le.0.) then
         fac=rdfac*1./rdfac
      else
         fac=rdfac
      endif

!c     first guess for ci: to determine wc,wj
      cint1=0.7*ca
      wc=vcm*cint1/(cint1+kc*(1.+o2/ko))
!c     wo=wc*o2*1000./(tau*cint1)
!c     wj=(wc-.5*wo)*pm/(wc+wo)/(1-gamma/cint1)
      wj=pm*cint1/(cint1+2*gamma)


      do flag=0,10
!c------- 1st case (left curve) ---- minwc ---------------
         if (wc.le.wj) then
            aaa=vcm
            bbb=kc*(1.+o2/ko)
            ddd=gamma
            eee=1.

            call psncalc(aaa,bbb,ddd,eee,gfac,gmin,ga,ca,  &
              humid,rd,fac,pnet,cint,gpres)

            if ((cint.gt.(cint1+.1)).or.(cint.lt.(cint1-.1))) then
               cint1=cint
               wc=vcm*cint1/(cint1+kc*(1.+o2/ko))
               wj=pm*cint1/(cint1+2*gamma)
            else
               goto 567
            endif

!c------- 2nd case (right curve) ---- minwj --------------
         else
            aaa=4*pm
            bbb=8*gamma
            ddd=gamma
            eee=4.

            call psncalc(aaa,bbb,ddd,eee,gfac,gmin,ga,ca,  &
               humid,rd,fac,pnet,cint,gpres)

            if ((cint.gt.(cint1+.1)).or.(cint.lt.(cint1-.1))) then
               cint1=cint
               wc=vcm*cint1/(cint1+kc*(1.+o2/ko))
               wj=pm*cint1/(cint1+2.*gamma)
            else
               goto 567
            endif
         endif
      enddo
 567  continue

      ci(j)=cint

      call calcTotals(j,pnet,gpres,rd,fac)

      return
      end   !subroutine psn6 end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine psncalc(aaa,bbb,ddd,eee,  &
          ggfac,ggmin,ga,ca,rh,rd,fac,pnet,cint,gpres)
!c
!c     analytical leaf photosynthesis simulation
!c     BALDOCCHI (1994. Tree Physiology 14: 1069-1079)
!c     validated by BALDOCCHI for:
!c           PAR:                               0-2000 [umol m-2 s-1]
!c           CA:                               50-650  [ppm]
!c           Tleaf:                             5-35   [deg.C]
!c           Boundary Layer Resistance (1/GA):  1-2000 [s m-1]
!c-----------------------------------------------------------------------

      USE constants
      implicit none

      double precision ca,ga,gpres,pnet,rd,cint,cs,rh
      double precision aaa,bbb,ddd,eee,ggfac,ggmin,fac
      double precision rrr,sss,ttt,QQQQ,PPPP,alfa,betta,gama,tetta,   &
          uuu,vvv,DDDD,pnet1,pnet2,pnet3,aquad1,bquad1,cquad1

!c     qqq       [mmol m-2 s-1]
!c     rrr,PPPP  [(mmol m-2 s-1)**2]
!c     sss,QQQQ  [(mmol m-2 s-1)**3]

      alfa=1600.*1600.*(1+ggmin/ga-ggfac*rh/1.6)
      if (alfa.ne.0.) then
         betta=1600.*ca*(ggfac*rh/1.6*ga-2*ggmin-ga)-   &
           1600.*1600.*(fac*rd*ggfac*rh/1.6)
         gama=ga*ca*ggmin+1600.*ga*fac*rd*ggfac*rh/1.6
         tetta=1600.*(ggfac*rh/1.6*ga-ggmin)

         rrr=(eee*betta+bbb*tetta-aaa*alfa+eee*alfa*fac*rd)/   &
           (eee*alfa)
         sss=(eee*gama*ca+bbb*gama-aaa*betta+aaa*ddd*tetta+   &
           eee*fac*rd*betta+bbb*fac*rd*tetta)/(eee*alfa)
         ttt=(-aaa*gama*ca+aaa*ddd*gama+eee*fac*rd*gama*ca+bbb*fac*rd    &
           *gama)/(eee*alfa)

         PPPP=(3*sss-rrr*rrr)/3.
         QQQQ=2.*rrr*rrr*rrr/27.-rrr*sss/3.+ttt
         DDDD=(PPPP/3.)*(PPPP/3.)*(PPPP/3.)+(QQQQ/2.)*(QQQQ/2.)
         if (QQQQ.lt.0) then
            vvv=-sqrt(abs(PPPP)/3.)
         else
            vvv=sqrt(abs(PPPP)/3.)
         endif
         uuu=QQQQ/2./(vvv*vvv*vvv)
         if (PPPP.gt.0) then
!c           eine reelle Loesung, zwei konjugierte
            pnet=-2*vvv*sinh(log(uuu+sqrt(uuu*uuu+1))/3.)-rrr/3.
         else
            if (DDDD.gt.0) then
!c              eine reelle Loesung, zwei konjugierte
               pnet=-2*vvv*cosh(log(uuu-sqrt(uuu*uuu-1))/3.)-rrr/3.
            else
!c              drei reelle Loesungen
               pnet1=-2*vvv*cos(acos(uuu)/3.)-rrr/3.
               pnet2=-2*vvv*cos(acos(uuu)/3.+2*pi/3.)-rrr/3.
               pnet3=-2*vvv*cos(acos(uuu)/3.+4*pi/3.)-rrr/3.

               if (alfa.gt.0.) then
                  pnet=pnet2
                  if ((pnet1.gt.0.).and.(pnet2.gt.0.).and.   &
                    (pnet3.gt.0.))then
                     pnet=min(pnet1,pnet2,pnet3)
                  elseif ((pnet1*pnet2*pnet3.lt.0.)  &
                    .and.(max(pnet1,pnet2,pnet3).gt.0.)) then
                     pnet=min(pnet1,pnet2,pnet3)
                  else
                     pnet=max(pnet1,pnet2,pnet3)
                  endif
               else
                  pnet=pnet3
               endif
            endif
         endif

!c        pnet  [umol CO2 m-2 s-1]
!c        cs    [ppm]
!c        gpres [mmol H2O m-2 s-1]
!c        cint  [ppm]

         cs=ca-pnet*1600./ga

         gpres=ggfac*1000.*(pnet+fac*rd)*rh/cs+ggmin
         cint=cs-pnet*1600./gpres
!c        print*,'kubisch: pnet cint gpres',pnet,cint,gpres
      else
!c        alfa=0: quadratische Loesung
         betta=-1600.*ca*ggmin-1600.*1600.*(fac*rd*ggfac*rh/1.6)
         gama=ga*ca*ggmin+1600*ga*fac*rd*ggfac*rh/1.6
         tetta=1600.*ga

         aquad1=eee*betta+bbb*tetta
         bquad1=eee*gama*ca-aaa*betta+eee*fac*rd*betta+aaa*ddd*tetta+   &
           bbb*fac*rd*tetta+bbb*gama
         cquad1=-aaa*gama*ca+eee*fac*rd*gama*ca+bbb*fac*rd*gama+aaa*   &
           ddd*gama

         pnet=real((-bquad1-sqrt(bquad1*bquad1-4.*aquad1*cquad1))/   &
           (2.*aquad1))
         cs=ca-pnet*1600./ga
         gpres=ggfac*1000.*(pnet+fac*rd)*rh/cs+ggmin
         cint=cs-pnet*1600./gpres
!c        print*,'quadratisch: pnet cint gpres',pnet,cint,gpres
      endif

      return
      end





!gfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine calcTotals(j,pnet,gpres,rd,fac)
!c
!c     init vphoto-common
!c-------------------------------------------===-------------------------
      implicit none

      common/ctrl1/needles,p_flag
      common/ctrl2/sigma,fac1
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)

      integer j,needles,p_flag
      double precision pnet,gpres,rd,sigma,fac1,ci,ga,pn,gs,gl,wd,fac

!c     fuer Koniferen muessen Gaswechselraten auf proj. intact twig area
!c     bezogen sein: umol/mmol pro m2 proj. intact twig area und s
!c     (also unter Beruecksichtigung von Zusammenballen der Nadeln
!c     in "intakte Zweige")
      if (needles.eq.1) then
!c        falls PSParam. von Raten pro proj.Flaeche bestimmt (p_flag=0)
         if (p_flag.eq.0) then
            pn(j)=pnet/2.57/fac1
            gs(j)=gpres/2.57/fac1
            wd(j)=rd/2.57/fac1
!        falls PSParam. von Raten pro Oberflaeche bestimmt (p_flag=1)
!        zB Picea
         elseif (p_flag.eq.1) then
            pn(j)=pnet/fac1
            gs(j)=gpres/fac1
            wd(j)=rd/fac1
         endif

!     fuer Blaetter muessen Gaswechselraten auf proj. leaf area
!     bezogen werden: umol/mmol pro m2 proj. leaf area und s, da
!     leaf area "projizierte leaf area" ist
      elseif (needles.eq.0) then
!        falls PSParam. von Raten pro proj.Flaeche bestimmt (p_flag=0)
         if (p_flag.eq.0) then
            pn(j)=pnet
            gs(j)=gpres
            wd(j)=rd*fac
!        falls PSParam. von Raten pro Oberflaeche bestimmt (p_flag=1)
         elseif (p_flag.eq.1) then
!ef         Spezialfall: Pinus penumbra Hyytiala, Loobos,
!ef         leaf area bleibt LAI, gas exch. per LSA wird auf LAI bezogen,
!ef         wd ist tagsueber 0.5*rdark, nachts 1.0*rdark
            pn(j)=pnet*2.7
            gs(j)=gpres*2.7
            wd(j)=fac*rd*2.7
         endif
      endif

      return
      end






!cgfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      function resp(t,d,e)
!c
!c     This function calculates the respiration by means of the
!c     ARRHENIUS-Equation
!c     T: Temperature of Leaf or Stem or Soil
!c     D: Scaling factor
!c     E: Activation Energy for Respiration
!c-----------------------------------------------------------------------

      USE constants

      implicit none

      double precision t,d,e,tk,rtk,resp

      tk=t+273.16
      rtk=rgas*tk

!c     if the parameters E and D = 0 then no respiration calculation needed
      if ((d.eq.0.).and.(e.eq.0.)) then
         resp=0.0
         return
      endif

!c     Arrhenius Equ.
      resp=exp(d-e/rtk)

      return
      end







!csoil--------------------------------------------------------------------
!c     PROXEL77 for TOPGAS            5.0.00              gy     21.09.00
!c------------------------------------------------------------------------
!c     Programmers PROXEL: Mauser, Joss, Reichstein
!c------------------------------------------------------------------------
!cinc~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!cinc  includefile soil1pix_dble.inc:
!cinc~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!cinc    integer nsl,nos
!cinc
!cinc    parameter (nos=1)         ! number of soils
!cinc    parameter (nsl=4)         ! number of soillayers + 1 !!!
!cinc
!cinc    integer nOfLayers,lowbouZero
!cinc    double precision pond,albedo,infill,
!cinc   +  runoff,actEvap,effPsi,interc,drainage
!cinc
!cinc    common/soiltyp/nOfLayers,lowbouZero,pond,albedo,infill,
!cinc   +  runoff,actEvap,effPsi,interc,drainage
!cinc
!cinc
!cinc    double precision  thick(nsl),zz(nsl),lambTemp(nsl),spHeatCap(nsl),
!cinc   +  SOM(nsl),Temp(nsl),hFlux(nsl),lambA(nsl),lambB(nsl),      ! som: read, but never used
!cinc   +  lambC(nsl),lambD(nsl),lambE(nsl),kSat(nsl),thetaR(nsl),
!cinc   +  thetaS(nsl),thetaMax(nsl,nos),vgN(nsl),vgM(nsl),vgAlpha(nsl),
!cinc   +  theta(nsl),psi(nsl),qH2O(nsl),
!cinc   +  qSatFlow(nsl),soilCond(nsl),cOrg(nsl),k1(nsl),k2(nsl),
!cinc   +  fractK1(nsl),respTP1(nsl),respTP2(nsl),respTP3(nsl),bd(nsl),
!cinc   +  respFP1(nsl),respFP2(nsl),respFP3(nsl),envFac(nsl),resp(nsl),
!cinc   +  rootDens(nsl),rootUp(nsl),psiRoot(nsl),bb(nsl)
!cinc
!cinc    common/soillayer/thick,zz,lambTemp,spHeatCap,bd,SOM,Temp,hFlux,
!cinc   +  lambA,lambB,lambC,lambD,lambE,kSat,thetaR,thetaS,thetaMax,
!cinc   +  vgN,vgM,vgAlpha,theta,psi,qH2O,qSatFlow,soilCond,
!cinc   +  cOrg,k1,k2,fractK1,respTP1,respTP2,respTP3,respFP1,respFP2,
!cinc   +  respFP3,envFac,resp,rootDens,rootUp,psiRoot,bb
!cinc~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



